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Abstract 

Population rate or activity equations are the foundation of a common approach to modeling 
for neural networks. These equations provide mean field dynamics for the firing rate or activity 
of neurons within a network given some connectivity. The shortcoming of these equations is 
that they take into account only the average firing rate while leaving out higher order statistics 
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like correlations between firing. A stochastic theory of neural networks which includes statistics 
at all orders was recently formulated. We describe how this theory yields a systematic extension 
to population rate equations by introducing equations for correlations and appropriate coupling 
terms. Each level of the approximation yields closed equations, i.e. they depend only upon 
the mean and specific correlations of interest, without an ad hoc criterion for doing so. We 
show in an example of an all-to-all connected network how our system of generalized activity 
equations captures phenomena missed by the mean field rate equations alone. 

1 Introduction 

Modeling the brain is confounded by the fact that there are a very large number of neurons 
and the neurons are heterogeneous and individually complex. Given current analytical and 
computational capabilities, we can either study neuronal dynamics in some biophysical detail 
for a small or medium set of neurons or consider a large population of abstract simplified neural 
units. We then can only extrapolate to the desired regime of large numbers of biophysical 
neurons. In particular, there is a dichotomy between network models that incorporate Hodgkin- 
Huxley or integrate-and-fire spiking dynamics and models that only include the rate or activity 
of neural units. While the rate description has yielded valuable insights into many neural 
phenomena it cannot describe physiological phenomena thought to be important for neural 
processing such as synchronization, spike-time dependent plasticity or any correlated activity 
at the spike level. Likewise, it is difficult to analyze or simulate a large network of spiking 
neurons. Our goal is to derive an intermediate description of neural activity that is complex 
enough to account for spike level correlations yet simple enough to be amenable to analysis 
and numerical computation for large networks. 

Rate or activity equations have been a standard tool of computational and theoretical neu- 
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roscience, early important examples being the work of Wilson and Cowan, Cohen and Gross- 
berg, Amari and Hopfield [Wilson and Cowan, 1972 Wilson and Cowan, 1973, Amari, 1975 



Amari, 1977[ [Hopfield, 1984[ [Cohen and Grossberg, 1983] , Models of this type have been 



used to investigate pattern formation, visual hallucinations, content addressable memory 
and many other questions [Ermentrout and Cowan, 1979 Hopfield, 1984 Ermentrout, 1998 



Bressloff et al., 2002 Coombes, 2005| . Naturally, these equations are so called because they 



describe the evolution of a neural activity variable often ascribed to the firing rate or synap- 
tic drive of a population of interacting neurons [Ermentrout, 1998 Gerstner, 2000[ . These 
equations are considered to represent the neural dynamics averaged over time or population 
of a more complicated underlying process. In general, these activity equations make an im- 
plicit assumption that correlated firing is unimportant. They are a "mean field theory" which 
capture the dynamics of the mean firing rate or activity that is independent of the influence 
of correlations, which in some cases may alter the dynamics considerably. As an example, 
the effects of synchrony, which have been proposed to be important for neural processing 
[Gray and Singer, 1989 Beshel et al., 2007] are not included. Here, we give a systematic pre- 
scription to extend rate models to account for these effects. 

An analogy for our problem and approach can be made to the field of equilibrium statistical 
mechanics. The statistics of such systems (e.g. the Ising model) in thermal equilibrium are 
described by a partition function, which is an integration over all configurations available to 
the system. For the Ising model this refers to all possible configurations of the individual spins. 
The partition function is akin to the generating function for a statistical distribution from which 
the moments or cumulants can be obtained. For the Ising model the first moment corresponds 
to the mean magnetization and the second moment describes the mean correlation between 
the spins. The linear response of the system is the magnetic susceptibility, which describes 
the reaction of the system to an external input. In general the partition function cannot be 
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summed or integrated explicitly. However, these moments can be obtained perturbatively 
by using the method of steepest descents to approximate the partition function. This then 
yields a systematic expansion and the lowest order is called mean field theory, since all higher 
cumulants are zero. By computing the expansion to higher order, the effects of correlations 
and fluctuations can be included. 

This procedure requires full knowledge of the underlying microscopic theory that is to be 
averaged over. In neuroscience, the underlying model is not completely known; it would require 
full knowledge of the different types of neurons, their membrane and synaptic kinetics, and 
their synaptic connectivity. However, given a particular mean field theory, one can ask about 
the minimal constraints this theory places upon the microscopic theory and its asymptotic 
expansion. Thus, although the full microscopic theory cannot be reconstructed, by constraining 
the expansion, the mean field theory can dictate the minimal structure of any extension of a 
set of rate equations. In this paper, we consider a well known neural rate equation and deduce 
the minimal structure we expect for a consistent extension which includes correlations. 

Buice and Cowan [Buice and Cowan, 2007| previously adapted a path integral formalism 
used in nonequilibrium statistical mechanics |Doi, 1976a[ Doi, 1976b, Peliti, 1985 to analyze 



the dynamics of a Markov model for neural firing. They derived a generating functional 
(expressed as an infinite dimensional path integral), which is specified by an "action" for the 
complete dynamical distribution of the model and showed that the mean field theory for that 
system corresponded to a Wilson-Cowan-type rate equation. They then analyzed the scaling 
properties for the correlations near criticality. They also showed how mean field theory could 
be corrected by using steepest descents to generate a systematic expansion that describes the 
effect of correlations. Here, we show that by taking explicit averages of the Markov model a 
moment hierarchy can be constructed. Each equation in the moment hierarchy is coupled to 
higher moments in the hierarchy. The hierarchy can be made useful as a calculational tool 



for the statistics of the dynamics if it can be truncated. We show that the moment hierarchy 
and the generating functional are equivalent and that the equations of the hierarchy are the 
"equations of motion" of the action in the generating functional. The truncation condition 
for the perturbation series of the path integral is also a truncation condition for the hierarchy. 
This provides for both a compact description of network statistics and a natural truncation or 
closure condition for a moment hierarchy. We can also show using the path integral formalism 
that the Markov model is a natural minimal extension of the Wilson-Cowan rate equation. 

Approaches to neural network modeling using statistical mechanics are not new [Hopfield, 1982 



Hopfield, 1984, |Peretto, 1984 Amit et al., 1985 . Those works were largely concerned with 
models adhering to detailed balance, whereas we make the explicit assumption that neural 
dynamics admits an absorbing state that violates detailed balance. In the absence of in- 
ternal activity and external stimulation, there will be no activity in the network. Other 
studies using a stochastic description of neural dynamics have considered the neurons in a 
background of Poisson activity with disorder in the connectivity [Amit and Brunei, 1997a 



Amit and Brunei, 1997b] , or considered neural activity as a renewal process [Gerstner, 1995 



Gerstner, 2000, Gerstner and Kistler, 2002 . Van Vreeswijk and Sompolinsky [1996,1998] demon- 
strated that disorder in network activity can arise purely as a result of disorder in the con- 
nectivity, without stochastic input. Kinetic theory and density approaches are investigated 
in [Nykamp and Tranchina, 2000 Cai et al., 2004, Ly and Tranchina, 2007| and mean field 
density approaches to the asynchronous state appear in [Abbott and van Vreeswijk, 1993 



Treves, 1993] . [Golomb and Hansel, 2000] study synchrony in sparse networks via a reduction 
to a phase model. Fokker-Planck approaches for networks appear in |Fusi and Mattia, 1999 



Brunei and Hakim, 1999 , N., 2000] [Brunei and Hansel, 2006 . Responses of single neurons 

driven by noise appear in [Plesser and Gerstner, 2000 Salinas and Sejnowski, 2002, Fourcaud and Brunei, 2002 

Soula et al., 2006] . Approaches to correlated neural activity including finite size effects appear 



in Ginzburg and Sompolinsky, 1994 , Mattia and Del Giudice, 2002 , Soula and Chow, 2007 , El Boustani and De 



In El Boustani and Destexhe, 2009 , the authors develop a moment hierarchy for a Markov 
model of asynchronous irregular states of neural networks which is truncated through a combi- 
nation of finite size and a scaling condition. Our work extends the results of [Ginzburg and Sompolinsky, 1994] 
by providing the systematic higher order expansion without explicitly requiring the considera- 
tion of the rest of the hierarchy. We also provide conditions for the truncation of the expansion 
and consider the network response to correlated input. Our expansion is not a finite size ex- 
pansion, although it can reduce to a finite size expansion under certain conditions (such as 
normalized all-to-all connectivity in the network). 

In section [21 we revisit the original Wilson-Cowan framework and propose a Markov model 
that has the minimal stochastic dynamics to produce the Wilson-Cowan equations. This 
will be more rigorously justified in section [H Section [3l presents the derivation of a moment 
hierarchy for this Markov model. After truncating, we provide a posteriori justification for 
the truncation. It will be seen in section [H that the validity of this truncation was in fact 
natural and did not require ad hoc assumptions. The truncation conditions turn out to be 
related to the proximity to a bifurcation point as well as the extent of connectivity in the 
network. We also make more precise the sense in which our Markov model is "minimal" 
by introducing the path integral formulation. The field theory formalism which appears in 



this paper arose in the context of reaction-diffusion problems. See [Janssen and Tauber, 2005 



Tauber et al., 2005[ for reviews of this formalism applied to reaction-diffusion and percolation 



processes. We demonstrate a simple example all-to-all system in Section [5] and show some 
simulation results. 
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2 Rate equations reconsidered 

We consider a population rate equation of the form 

dta{x,t) = -aa{x,t) + f (^j d'^y w{x,y)a{y,t) + I{x,t)^ (1) 

where a{x, t) is a measure of local neural "activity" at location x G W^, q is a decay constant 
(often equated with the membrane time constant or a synaptic time constant), /(s) is the firing 
rate or gain function describing how input affects the activity, I{x, t) is a time dependent 
external input to location x, and w{x,y) is a weight function describing how a neuron at 
location y affects neurons at location x. Equation ([1]) is the standard form of rate equation seen 
in the Wilson-Cowan equations [Wilson and Cowan, 1972 Wilson and Cowan, 1973] . While 



we use this form in our paper, our results can be adapted to any other type of rate or activity 



equation. The exact nature of the activity a{x, t) is open to interpretation [Gerstner, 1995 



Ermentrout, 199"8] . It could be envisioned as the time average, population average or ensemble 



average of neural firing or synaptic activity. In any case, the picture is that of some underlying 
process whose degrees of freedom have been marginalized to generate an effective theory with 
simpler variables. 

We imagine that the typical rate equation is produced by some marginalization process over 
both disorder and extra degrees of freedom. Hence, it may be possible to derive a generating 
function for the statistics of the marginalized process. The lowest order in the steepest descent 
expansion of the generating function describes "mean field theory" , which gives the rate equa- 
tion. Since the operation of marginalization is dissipative, we cannot recreate the underlying 
microscopic process exactly with only the rate equation alone. However, the mean field theory 
places constraints upon the structure of the dynamics, enabling us to investigate the structure 
of higher order statistics implied by the structure of mean field theory. In the original deriva- 
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tion by Wilson and Cowan Wilson and Cowan, 1972 , Wilson and Cowan, 1973 , the activity 
variable was presumed to describe the fraction of neurons firing per unit time within some re- 
gion of the brain. There are two main features of this interpretation which bear emphasizing. 
First, the rate equations were originally understood to be equations providing the dynamics of 
the probability that a neuron at x will fire at time t. There is therefore an implied underlying 
probabilistic model. Second, the probability a{x, t) applies to all neurons within some region 
of the brain, not just a single neuron. Thus, there is a spatial averaging component implicit 
in the equations. The original Wilson-Cowan rate equations thus described the dynamics of 
the probability for a neural aggregate in the brain. Another feature implicit in the Wilson- 
Cowan equations is that these probabilities are independent for each neuron. This implies that 
the Wilson-Cowan picture is one in which neurons fire with Poisson statistics with firing rate 
determined by a{x,t), a picture supported by neural recordings [Softky and Koch, 1993| . 

Given this perspective, one might consider what processes may underly rate equations. 
One route is to treat the fundamental, small-scale dynamics as a probabilistic process, for ex- 
ample a Markov process. In this case, the basic description for neural activity will be provided 
by a master equation governing the evolution of probabilities for different neural configura- 
tions. This route obscures the source of uncertainty in neural activity in favor of directly 
modeling the probabilistic activity. This tactic has been used to model the so-called asyn- 
chronous irregular states seen in some neural models [Van Vreeswijk and Sompolinsky, 1996 



El Boustani and Destexhe, 2009| . Another route would be to employ the strategy of kinetic 



theory [Nicholson, 1992 Ichimaru, 1973| and define a continuity (i.e. Klimontovich) equation 
for the probability density of a network of deterministic neurons. For an example of this ap- 
proach applied to coupled oscillators, see Hildebrand et al., 2006 Buice and Chow, 2007| . In 
that work, the probabilistic aspects of the model arise from the distribution of driving frequen- 
cies and initial conditions. Ultimately, the difference in the two approaches is the origin of 
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stochasticity, i.e. whether it is imphcit in the dynamics of the neurons or an emergent property 
of the interaction of deterministic neurons (e.g. chaos). In either case, the final product is an 
effective stochastic dynamical system. In this paper, we will follow the approach of assuming 
an underlying probabilistic model given by a master equation, so that any emergent chaos 
has already been absorbed into the dynamics. We then seek a minimal stochastic model that 
will produce the Wilson-Cowan rate equation at the mean field level. We can then formulate 
equations governing the fluctuations of this model. In this section we motivate such a minimal 
model qualitatively, leaving a more rigorous approach for section HI 

Our primary interest is in tracking the statistics of active neurons. A simple master equa- 
tion whose mean field statistics for neural activity is represented by the Wilson-Cowan rate 
equations is given by 

^^^^1^ = '^[a{ni + l)P{ni+,t) -aniP{n,t) 

i 

+ F, (fii.) P{n,_,t) - Fi (n) P(n, t)] (2) 

where -P(n, t) is the probability of the network having the configuration described by n = 
{ni,n2, • • •} at time t, and rij is the number of active neurons at location i. Neurons relax 
back to the inactive or quiescent state with rate a, which appears as a decaying transition 
in the master equation. Configurations nj+ and denote the configuration n where the 
ith component is fij it 1, respectively. The rate at which a neuron at location i becomes 
active is given by the firing rate or gain function Fi{n)^ which is an implicit function of 
the weight function Wij and external inputs /j. One of the crucial elements of the ensuing 
calculation is making the connection between the gain function Fi{n), which appears in ^ 
and f{Ii{t) + Ylj WijUj), which appears in ([1]). 

In general, we cannot solve for t) in ([2|) explicitly. One strategy is to derive an 
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expansion of P{n,t) in terms of its moments {ni{t)nj{t')nk{t") ■ ■ ■) where the expectation 
value is over all statistical realizations of the Markov process. The first moment 



(t) = {n,{t)) = ^niP{n,t) 



(3) 



n 



is a measure of the mean activity in the network. We obtain an equation for ai{t) by multiplying 
equation ([2]) by rij and taking the sum over all configurations n. However, this equation is 
not closed (i.e. it depends upon the second and possibly higher moments). An equation for 
the second moment can be similarly constructed by multiplying ([2|) by riirij and summing 
over all configurations. The resulting equation will depend on the third and higher moments. 
Continuing this process will result in a moment hierarchy with as many equations as there 
are locations, which could be infinite. In general no finite subset of this hierarchy is closed. 
This means that if we wish to have a closed set of equations then we need to make some 
approximation which allows us to truncate the hierarchy. 

The simplest way to close or truncate the moment hierarchy is to assume that all the 
higher order moments factorize into products of ai{t). This is the naive mean field as- 
sumption where all cumulants are zero. For example, the second cumulant (i.e. variance) 
{ni{t)nj[t')) — ai{t)aj{t') is set to zero. For our master equation ([2]), this assumption yields 
(see the computation in the next section) 



which is similar in form to the rate equation ([T]) for a discrete domain and with the firing 
rate function given by Fi{a). However, the fact that statistics in the brain are observed to 
be near Poisson is devastating for this naive mean field assumption because every cumulant 




(4) 
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is comparable to the mean, implying any such truncation of the resulting hierarchy is not 
justifiable. 

Here, we describe an alternative means of truncating the hierarchy consistent with near 
Poisson firing statistics. In this case, we observe that the equation for ai{t) can be written as 



where g is some functional dependent upon the higher moments in the hierarchy. In order to 
choose a reasonable approximation and get a finite system of closed equations, we must identify 
a finite set of higher moments in terms of which we can express the remaining moments. We 
are guided by the indication that neuron firing statistics are near Poisson. Not coincidentally, 
the solution to the master equation in the case where the gain function Fi is constant or linear 
is exactly Poisson with mean rate (i.e. stochastic intensity) determined by the Wilson-Cowan 
rate equation. In order to truncate the hierarchy, we perform a change of variables to measure 
the deviations of each cumulant from the value under Poisson statistics. This new hierarchy is 
truncatable, as will be demonstrated a posteriori. From the perspective of solving the master 
equation, this new hierarchy is the natural one because the underlying statistical model is a 
point process. 

The moment hierarchy approach does not make any approximation. It is a change of 
variables from the distribution P{fi, t) to moments of that distribution. The approximation 
arises when we truncate this hierarchy in order to render the equations tractable. The simplest 
truncation is mean field theory. The first order corrections to mean field theory are given by 
truncating at the next order. Truncation of the moment hierarchy requires some justification. 
We will demonstrate below that this justification in the neural case may be provided by the 



dtai{t) 



aai{t) + Fi (a) 



+9 [{'rii{t)nj{t)), {ni{t)nj{t)nk{t)), {ni{t)nj{t)nk{t)ni{t)), ■ • •] 



(5) 
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large spatial extent of neural connectivity and the distance of the system from a bifurcation. 



3 Truncation of the Moment Hierarchy 

We will derive a moment hierarchy from our master equation ([2]) and then show how it can be 
truncated. To get an equation for the first moment ai{t), we multiply the master equation ([2]) 
by rii and sum over all configurations n: 



^^^dP{^^ ^ I «(^fc + l)P{nk+,t) - ankPjn, t) 

n ft K k 

+ Fi (ni_) Pini.,t) - Fi (n) t)} (6) 



The first two terms on the right hand side simplify to 

I + l)P(nfc+,t) - anfcP(n,t) I = a^'^Ui {nkP{n,t) - nkP{n,t)) 

ft \ k / k^i ft 

+a^{ {ui - 1) niP{n, t) - njP{n, t) ] 
ft 

= —a''^^niP{n,t) = —aai{t) (7) 



n 



The first equality results from re-indexing the summation over from (0, oo) to (l,oo). We 
leave the summation indicated as over all configurations n because the factor of prevents 
the term from contributing. We have also separated out the terms where i = k; the only 
term which survives is one of these. Note that we have made no approximations thus far. The 
terms involving the function Fi{n) take the form 



^ ^ {Fi (Hi.) P{ni-,t) - Fi (n) P{n, t)} 
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+ ^ (n, + 1) Fi{n)P{n, t)-J2 niF,{n)P{n, t) 



Y,m)P{n,t) (8) 



Unlike the first term, we cannot directly represent this term as a function only of ai(t), due to 
its nonlinear nature. The equation for the mean is therefore: 

^ = -aai + (F, (n)) (9) 

where {Fi{n)) will include higher order moments. 

We continue by constructing an equation for the second moment 

Nij = {niUj) = Y^ninjP{n,t) (10) 
n 

and third moment 

Nijk = {riinjUk) = Y^ninjnkP{n,t) (11) 

n 

to obtain 

^ = -2aNij+a5ijai + {njFi + niFj) + 5ij{Fi) (12) 
= -3aNijk + a{ninj6irn + ninrn6jrn + njnrnSij) - a6ij6jrn{ni) (13) 
+ {niUjFk + njnfcFj + UjUkFi) + {njFkSik + UkFiSij + UiFjdkj) (14) 

3 (tJ, j i<i ) (Jjm, 

Since we expect solutions to be near Poisson, we transform the hierarchy to describe the 
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departure of moments from Poisson statistics. For a Poisson distribution, the cumulants are 
all equal to the mean of the distribution. Hence, we introduce what are called "normal ordered 
cumulants" , which measure the "deviations" of the cumulants from Poisson values. The first 
normal ordered cumulant is the same as the first moment ai{t). The next two are given by 

Cij = Nij - aiUj - aiSij (15) 

and 

Cijk = Nijk — Nijttk — Njktti — Nikttj + 2aiajak 

-{Nij - aiaj)5jk - {Nik - aiak)Sij - {Njk - ajak)6ik + 2ai6ij6jk (16) 

The normal ordered cumulants can be computed using a recursive algorithm. The algorithm 
involves replacing all moments Nijk... with "normal-order-corrected" moments recursively by 
subtracting terms with coincident or "contracted" indices which reduce the order of the mo- 
ment. For example, the ordinary second cumulant is simply iVjj — OjOj. To compute the normal 
ordered version, we replace the moments appearing in the expression with the normal-ordered- 
corrected forms, i.e. set Nij — > Nij — aiSij. The term subtracted results from the contraction 
of the i,j indices, i.e. Nij Ni5ij = aiSij. The third cumulant is more complicated but 
follows the same strategy. The important thing to note here is that the higher moments must 
be independently corrected for each group of contracted indices. The ordinary third cumulant 
is given by 

{{rii - ai) {rij - aj) {uk - ak)) = Nijk - NijUk - NjkUi - Nikaj + laiUjUk (17) 
To obtain the normal ordered cumulant, we first make the replacement Nijk N^jk — Nijdjk — 
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^ik^ij — Njk^ik — CLi6ij5jk- We then must correct for all appearances of the second moment 
Nijn resulting in (jl6p . This algorithm systematically removes the underlying Poisson contri- 
butions (at all tensor ranks) and leaves us with the normal ordered cumulants. For a Poisson 
distribution, all Cijk— = except the first, aj. An alternative rationale for normal ordering 
will appear in Sec. [H 

Transforming the first three equations of the hierarchy ()12p . and (|13p yields 

= -aai + {F,) (18) 

dC- ■ 

= -2aCij + {{rij - aj)F,) + {{m - a,)Fj) (19) 

dC- -fc 

= —SaCijk + {{ui — ai){nj — aj)Fk + permutations) (20) 

where we must re-express the expectation values involving Fi in terms of the normal ordered 
cumulants. 

Since Fi{n) is defined in terms of the vector of active neuron numbers, n, its expectation 
value will be naturally expressed in terms of the moments, Nij^..., as given by the Taylor 
expansion of F[n). 

mn)) = E + ^ E PtN,k + ■■■ (21) 

We have implicitly defined the notation that F^^ is the derivative of Fi with respect to 
nj,nk,-''- Note that this expansion only applies to the expectation value of Fi{n). We 
need to re-express this series as an expansion in terms of the normal ordered cumulants. 
This transformation of variables will rearrange the terms and result in a new series with 
new coefficients that sums to the same result as the original series. For example every term 
proportional to Nij^... will contribute a term proportional to a^ajak • • • to the normal ordered 
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cumulant expansion. This means that we can write the expansion in the form 



(Fi) = Fi{a) + ■ ■ ■ 



(22) 



as we would expect. However, there are also contributions from the normal ordering correc- 
tions. The simplest are those which arise due to every index being coincident or contracted at 
each order. This produces a term linear in a^. These corrections form the series 



where by F^"" we mean the mth derivative of F{n) with respect to nj. Were -Fj(n) a polynomial 
this procedure would truncate at the highest order of the function. However, for an arbitrary 
general function these corrections quickly become unwieldy as one proceeds through the orders 
of the expansion to include corrections to the terms that go as aiaj, aiUjak ■ ■ ■■ 

Our perspective has been to interpret the Wilson-Cowan equation as describing the mean 
field of some Poisson process with activating and decaying transitions. Hence, the Wilson- 
Cowan equation should be the mean field solution of the normal-ordered cumulant hierarchy 
(|18p . However, the gain function in the Markov equation is not the same as the gain function 
in the Wilson-Cowan equation. The Wilson-Cowan gain function is the normal ordered version 
of the Markov gain function. Thus, in order for the Wilson-Cowan gain function to have the 
form of f{si), where Sj = Y^- Wijaj + we assume that this re-summation produces 




(23) 



m=2 j 



{F,) = f{si) + h.o.t. 



(24) 



where and the higher order terms (h.o.t.) are dependent upon the higher normal ordered 
cumulants according to the Taylor series expansion of f{s). It will be seen in Sec.[l]that there 
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always exists a Master equation gain function F such that / is expressible in this form and 
the resummation works to produce the same /(s) (and derivatives thereof) at every order in 
this expansion. 

We now return to the series (j23|) to consider terms with precisely one factor of Cij{t). 
At each order m in the series (i.e. the order which contains the mth moment), there are 
m(m — l)/2 terms which have one factor of Cij{t). These terms sum to give 

(\ m—2 
Y^w.Mt)] Y.'^ikWuCM (25) 
i / kl 

which can be rewritten as 

\f" Yl ^i^t) + w^3^ikCjk{t) (26) 

\ i / jk 

Other terms are at least second order in Cjj, higher normal ordered cumulants, or corrections 
from normal ordering. Our first generalized activity equation excluding terms dependent on 
third and higher normal ordered cumulants is therefore 

dtai{t) = -aaiit) + f ^y^^Wijajjt) + Ii{t) 

Yl ^i^t) + Y Wij^ikCjkit) (27) 

\ J / jk 

We can take this same approach in order to compute an equation for Cij{t) and obtain 
j^Cij{t) = -2aCij{t) + /' {si) Y mkCkjit) + /' (sj) Y WjkCkiit) 

k k 

+/' [si) Wijajit) + /' [sj) WjiUiit) (28) 
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Equations (j27p and ()28p constitute a closed set of equations for the mean and variance of a 
Wilson-Cowan network of neurons. ai{t) represents the Poisson rate and Cij{t) measures the 
deviation of the variance from Poisson statistics. These equations are the minimal consistent 
extension of the Wilson-Cowan rate equations to include the effects of higher order statistics. 
Higher order corrections can be incorporated by adding terms involving higher order cumulants 
into (j27p and (j28p and including equations for these higher order cumulants. We also note that 
the gain function need not be analytic everywhere for this expansion to work. It can contain 
a countable number of non-continuous or non-differentiable points. The equation would be 
corrected with the inclusion of impulse function terms at these singular points. 

An immediate noteworthy consequence is that Cij{t) will only have substantial input when 
the activity is such that f'{s) is large. As an example suppose f{s) is a simple sigmoid function. 
In this case, f'{s) is peaked at threshold (where we define threshold to be the central point 
of half maximum) and zero far away from threshold. Reasonably, we have the result that 
correlated activity will only increase when the input to a neuron is near threshold. If the slope 
of the sigmoid is such that f{s) is a step function, or near to a step function, then Cij{t) will 
receive input only when the activity is precisely near threshold. Also notice that the strength 
of the input to Cij{t) is proportional to the weight Wij between the neurons in question as well 
as the mean activity. An initial check on the equations is that Cij decouples from Oj in the 
case where /(s) is linear or constant. 

To consider the dynamics of large scale neural activity, we can take the continuum limit of 
these equations to get equations for the mean activity a{x,t) and correlation C{x,y,t). 
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X / dydzw{x,y)w{x, z)C{y, z,t) 



(29) 



and 



di 



Cix,y,t) = -2aC{x,y,t) + f {s{x)) J dzw{x,z)C{z,y,t) 
+/' is{y)) j dzw{y,z)C{z,x,t) 
+/' (•5(3^)) w{x, y)a{y, t) + /' {s{y)) w{y, x)a{x, t) 



(30) 



These are the generalized activity equations. Had we wished to include even higher moments 
we could have continued through the hierarchy. For simplicity of illustration, we truncate the 
hierarchy at this level. 

3.1 Criticality and truncation of the hierarchy 

Although we have derived equations for the mean activity and equal-time correlation, there 
are some outstanding issues. The primary concern which must be addressed is that we require 
some justification for the truncation of the hierarchy at the level of the two-point correlation 
function Cij{t) instead of allowing higher moments to interact with the mean activity. 

Consider the mean field equation without the correction due to correlated activity (jH). 
Define to be some steady state solution to this equation and linearize equation around 
this solution. The perturbations (5aj(t), from this steady state solution obey the equation 




(31) 
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We rewrite this equation as 

dt5a,{t) = -Y,^ij[a'^]6aj{t) (32) 
j 

where the matrix Tij is defined by 

r,, [aO] = a6i, - f ^ w^ja] + h (t) j (33) 

If all of the eigenvalues of are positive, then the solution a-* is stable. Likewise, negative 
eigenvalues indicate instability. Criticality is the condition of marginal stability, in which one 
or more of the eigenvalues are 0. 

Returning to the equation for Cjj(t), we see 

jCij{t) = -Y,{'^^kW)\Ckj{t)+Tjk[a{t)]Cik{t)) 

+/' {si) Wijaj{t) + /' (sj) Wjittiit) (34) 

We assume that the mean field solution is stable. In addition, we assume, per the trunca- 
tion hypothesis, that the steady state value of Cij{t) does not appreciably alter either Oj or, 
therefore, the matrix Tij. In this case, Tij has all positive eigenvalues and is diagonalizable. 
Define 

Kj = XiSij = UaTikU^l (35) 

Ik 

to be the diagonalization of r^. We also define the shorthand 

^ij = f {si) Wijaj{t) + /' (sj) Wjiai{t) (36) 
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for the driving terms in equation ([28]) . The steady state solution is given by 



■kl 



(37) 



llkk 



Notice that each term contributing to the magnitude of Cj'j is attenuated by a sum of eigen- 
values. The magnitude of the eigenvalues Aj determines the distance of the system from a 
bifurcation or criticality, i.e. it is the distance from the onset of an instability. Thus, the fur- 
ther the system is from criticality the more attenuated the fluctuations and the more justified 
we are in truncating the hierarchy. Conversely, the closer the system is to criticality the more 
the approximation breaks down. At criticality, this solution (j37p becomes singular. This is 
an indication that criticality is a fluctuation dominated, as opposed to mean field dominated, 
regime. A similar argument will extend to any equation in the hierarchy and we are left with 
an intuitively satisfying result: stability smooths out fluctuations. This argument is what 
allows us to disregard the effects of still higher moments upon the mean activity and truncate 
by considering only the two-point correlation function's effect on the mean. 

It is also worth noting that the eigenvalues relevant for the dynamics of the two point 
correlation Cij are the sums of the eigenvalues of the mean field equation, Aj + Xj. In the case 
that is stable, not only will Cij be stable but it will relax to equilibrium faster in general 
than the mean field solution. In kinetic theory, this is akin to the Bogoliubov approximation, in 
which the collision term is computed by solving for the steady state of the two point correlation 
on the assumption that the correlation function reaches steady state on a time scale shorter 
than the mean field. In our case, this approximation leads to 
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WijUj 



it)+Iiit)]^mjWikC% (38) 

jk 



Note immediately that the input from correlated activity decouples in the case that the firing 
rate function is in a linear or constant region, since the coupling is proportional to the second 
derivative of f{s). 

The extent of neural interconnections also has an effect on the size of correlations and the 
ability to truncate the hierarchy. Consider that neuron i has Ni pre-synaptic neurons (i.e. 
neurons for which Wij 7^ 0). Further, let the average connectivity weight over all inputs be 
wq. In this case we can approximate Wij ~ wo/Ni. The steady state value of Cij{t), which is 
determined by the driving term Aij{t), is seen to be inversely proportional to the number of 
pre-synaptic neurons due to the linear dependence on Wij of Aij . In the most extreme case the 
number of pre-synaptic neurons is the entire network, so Ni = N, and the correlation function 
Cij{t) becomes simply a finite size effect, going as Smaller system sizes in general will 

have larger correlations, which is intuitively sound. More generally, as long as we can bound 
the total input to any given neuron, we can define = minA'j and scale all weights so that 
they can be written 

N.m WM 

"^'^ IFW ^^^^ 

where wm is the maximum total input to any given neuron. This allows us to treat Nm as 
an effective system size. Larger N^a reduces the effects of fiuctuations at a given distance to a 
bifurcation. 

We have two competing effects. On one hand, we have the system size governing the 
magnitude of correlations. On the other, the distance of the system to a bifurcation likewise 
affects the size of fluctuations. The relative tradeoff of the two, from the definition of C?- 
is given by the product of the smallest eigenvalue of Tij and the number of pre-synaptic 
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neurons Ni. We will demonstrate this relationship more precisely when considering the all-to- 
all network in section [5j 

Finally, it is worth pointing out the effect of input upon the hierarchy. If the input is 
another Poisson process then the only equation which is affected is the equation for ai{t). 
The higher equations in the hierarchy are only affected by this input through its effect on 
the mean activity ai{t) and the firing rate function /(s). In general, this suggests that large 
external inputs will actually reduce fiuctuations, depending on the form of /(s), in the sense 
of driving the system towards Poisson-like behavior, a reasonable result. In particular, if /(s) 
is a saturating function, then the correlations will decouple from the equation for ai(t) and 
the source terms for higher correlations will be driven to zero, leaving the system described 
completely by the rate equations. The analogous situation for a ferromagnet is driving the 
system with a large external magnetic field. 



4 Path Integral Solution of the Master Equation 

We have thus far demonstrated how a minimal Markov model consistent with the Wilson- 
Cowan rate equation can be used to derive generalized equations in a hierarchy of moments. 
Although we truncated this hierarchy at second order, one can in principle truncate at any 
desired cumulant, although the calculations become successively more cumbersome. Here we 
show that the moment hierarchy is equivalent to a path integral or field theoretic approach, 
which systematizes the perturbation theory for the statistics of the network by providing rules 
for the construction and evaluation of the cumulants. Another major benefit is that it provides 
a systematic means for obtaining moment truncations or closures. The path integral represen- 
tation of the master equation ^ was derived by Buice and Cowan [Buice and Cowan, 2007| by 
modifying a method originally developed for reaction diffusion systems |Doi, 1976a Doi, 1976b 
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Peliti, 1985 . We quickly review the representation and then detail how the generalized equa- 
tions can be derived from this representation. 

The moment generating function for the probability density P{fi, t) is given by 



ft \ i / 



(40) 



where the sum is over all configurations of n. Moments of P are obtained by taking derivatives 
of the generating function with respect to Ji. For example {niUj) = d'^ Z/dJidJj\j^Q. The 
cumulants can be obtained by taking the derivatives of W[J] = \nZ. Field theory general- 
izes the generating function over a set of discrete variables to a generating functional over 
functions or fields. The result is a functional or path integral over all possible paths of time 
evolution for the system, weighted by the probability of that particular evolution. While it 
is sometimes possible and desirable to take the spatial continuum limit of the neural system, 
this is not necessary for the path integral approach. Here we use continuum spatial notation, 
although the results carry through in the case where x indexes a discrete variable. Buice 
and Cowan [Buice and Cowan, 2007| showed that the generating functional for the master 
equation ([2]) is given by 

Z[J{x, t), J{x, t)] = e^W^'*).-^>'*)] = J V^Vipe-^^^^'P^e^-^+'P-^ (41) 



where 



S[(p,ip\= I (fxdt 



d 



W[^{xM (42) 



is called the action and we use the notation u ■ v = f d'^xdtu{x,t)v{x,t). The fields ip and (p, 
which are obtained from the configuration variables n, are defined below [Janssen and Tauber, 2005 



Tauber et al., 2005| . The asterisk denotes convolution of the weight function with the inputs 
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and the term VF[(^(x,0)] is the cumulant generating functional of the initial distribution and 
takes into account arbitrary distributions in the initial condition. For example, if the initial 
state is described by Poisson statistics, we have VF[(^(x,0)] = / d'^x ao{x)ip{x,0), where ao(x) 
is the mean of the Poisson distribution at x. Analogous to the generating function for discrete 
variables, functional derivatives with respect to J(x, t) yield the normal ordered cumulants, 
such as 



6J{x,t) 



W[J{x,t),J{x,t)] 



{ipix,t)) = a{x,t) 



(43) 



J, ,7=0 



and 



6J{x,t) 6J{x',t') 



W[J{x,t),J{x,t)] 



J,J=0 



{^ix,tMx',t'))-{ip{x,t)){ip{x',t')) = C{x,t-x',t': 

(44) 



Within this formalism |Doi, 1976a Doi, 1976b Peliti, 1985 Buice and Cowan, 2007| , expecta- 
tion values of products of (p are the normal ordered cumulants found in the moment hierarchy. 
The normal ordered cumulant C{x,y,t) from (jlSp results from setting t = t' in C{x,t; x' ,t'). 
The field tp{x,t) is a "response" field and expectation over functions of it yield Green's func- 



tions or response functions for the dynamics [Martin et al., 1973 . The Ito convention is taken 
for the Langevin equation (j45p so that moments that involve combinations of (p and if that 
are evaluated at the same time are zero. More specifically, the convention taken is such that 
{(f{x,t')^p{x,t)) = Oiit<t'. 

We can heuristically derive the action ()42p , which was derived explicitly in [Buice and Cowan, 2007| 
and show that it represents a minimal model of the Wilson-Cowan equation where the activity 
is to be interpreted as a stochastic intensity or rate of a Poisson process. Consider an effective 
Wilson-Cowan Langevin equation 



dtn = —an + F{n) + no{x)6(t — Iq) + ^{x, t) 



(45) 
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where n{x,t) is the neural activity at location x and time t, £,{x,t) is an effective stochastic 
forcing with probability density functional and the firing rate function F is not necessarily 
that which appears in equation ([T|), but rather is that in the master equation We will show 
that Poisson noise is necessary to match the Buice and Cowan action (|42p . The probability 
density functional for n{x, t) can be written formally as 



P[n] oc j Vi5 



d 

—n + an — F (n) — nQ{x)5{t) — ^ 



P[i] (46) 



where 5[-\ is the functional generalization of the Dirac delta function. The probability density 
(j46p constrains the field n{x^ t) to obey the Langevin equation (j45p with initial condition no(a;). 
The delta functional is defined by the generalized Fourier transform 

PH c / nV^^^ (- / ,'.,tHa,n + o„ - F(„) -„„)-/ ^) PKI (47) 

and n is integrated along the imaginary axis. We can now integrate over the stochastic variable 
^ to obtain a noise generating functional defined by 

.«'I^1./P«»p(-/A..«)PK, (48) 

We choose ^ such that 

W[n] = (e" - 1 - n)F + an{e-^ - 1 + n) (49) 
which is consistent with a Poisson activation at rate F and a Poisson decay at rate a. We 
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next transform to the new variables 

if{x,t) = n(x,t)e-"(^'*) 

^{x,t) = e^(^'*)-l (50) 

The transformation (j50p to the new fields serves to simplify the noise generating functional 
and results in an action that has the form of (I42p but with a different gain function. This new 
action is reconciled with (j42p by the normal ordering operation. The reason this is necessary 
is because the Ito convention used to interpret the action would not hold uniformly for the 
transformed fields since in performing the transformation (jSOp . the <p and ip fields inside the 
gain function are evaluated at the same time and moments between these particular instances 
of the fields would not necessarily be zero. This inconsistency can be corrected by redefining 
(i.e. normal ordering) the terms in the gain function. As an example, consider the firing 
rate function to be F{n) = {w ■ nf' . After transforming, it becomes F = {w ■ {(pip + (/j))^, 
which mixes response and configuration operators at the same time point. To restore Ito 
convention, we normal order so that response and configuration variables are no longer mixed. 
We do this by considering the n(rE, t) operators at separate times t, t' and computing how 
the operators approach each other as t — > t'. The properties of the response field provide 
\\m.t^t'+{'P'{x-,t)(p{x\t')) = 5{x — x') and \\mt^t'-{'P'{x-,t)(p{x\t')) = so that F becomes 

lirn((zi; • {(p{x,t)(p{x,t) + (p){x,t)){w ■ {Lp{x,t')(p{x,t') + (p){x,t'))) = 

{{w ■ {^{x, t)ip{x, t) + ip){x, t)f + ■ {^{x, t)ip{x, t) + v?)(x, t))) (51) 

Hence, restoring the Ito convention requires the replacement n? v? + n (albeit in the 
variables 99, ip and similarly for higher powers of n). This normal ordering will adjust the gain 
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function to be f{w * {ifxp + ip)), leaving us with the action (j42p . From the perspective of the 
original master equation ([2]), the transformation to ip,(p is equivalent to expanding solutions 
to the master equation around a Poisson solution. 



4.1 Closed activity equations from the path integral 

We derive the generalized activity equations for the normal ordered cumulants directly from the 
generating functional by Legendre transforming to an effective action and then calculating the 
extrema of the effective action [Zinn- Justin, 2002 Cornwall et al., 1974, Buice and Cowan, 2007| . 



We first perform the computation for the mean activity a{x,t) and then show how to gener- 
alize to arbitrary numbers of cumulants. The generating functional (|41|) can be written more 
compactly as 



exp{W[J^'ix,t)]) 



j V^^ exp (^-5[$^] + j <fx dt Jf'ix, t)^ (52) 



where we define ^fj,{x,t), where fj, G {1,-1} such that <I>i(x,t) = (p{x,t) and $„i(x,t) = 
(f{x,t). We use Einstein summation convention (i.e. when the same index appears twice (one 
upper, one lower) in equations, summation will be implied). Similarly we define J^{x,t) via 
Ji = J and J-i = J. We can "raise" an index ^ via multiplication with the matrix 



V 



1 

1 



(53) 



so that = J and J ^ = J. Thus we have J^$^ = J^<I>i + J ^^-i = Jip + J(p. In order to 
streamline our notation, we will also define the dot product as the integral 

■a^ = j (fxdt J^(x, t)a^{x, t) (54) 
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For functions of more than one spatial variable, this inner product notation will generalize to 
the trace of the matrix product, i.e. 

A{x,y,t) ■ B{x,y,t) = J dtdxdyA{x,y,t)B{x,y,t) (55) 

We will also write the action as 5'[<I>^] = f d'^x L[<I>^], where L is the integrand of the action 
in dlS]). 

Define = (<J>^). The "effective action" r[a^] for is derived by a Legendre transforma- 
tion 

r[a/,(x, t)] = Jf'ix, t) ■ a^^ix, t) - T^[J^(x, t)] (56) 



where the conditions 



'IM^iil = (58) 

dafj,{x,t) 

are enforced. In analogy with classical mechanics, the extrema of the effective action 



5T[a^{x',t')] 
Sa^{x,t) 



(59) 



give the equations of motion or the activity equations for a^{x,t). 

In general, we will not be able to compute the equation of motion exactly since the path 
integral in (I52p cannot be computed exactly. However, we can perform a steepest descent 
asymptotic expansion of (j52p and compute the activity equation perturbatively. In field theory, 
this is known as the loop expansion because the terms in the expansion can be represented 
by Feynman diagrams whose order is given by the number of loops that diagram possesses. 
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Substituting for M^[J'^] using the Legendre transformation ([56]) gives 

exp(-r[a^]) = J exp (-5[$^] + • - • a^) (60) 

where we have suppressed the x and t arguments. Defining a new variable = — and 
using (j58|) gives 

exp(-r[a^]) = J V^^ exp {-S[^^ + a^] + • (61) 

where we define r^[a/i] = 5T/5a^. We now expand 5'[\I'/i + a^] in a functional Taylor series to 
obtain 

S[^^ + a^] = S[a^] + L^K] • + ^L^^'ia^] ■ + V[^,,, a^] (62) 

where represents the functional derivative of L[a^ with respect to and similarly for L^*^ 
and higher derivatives. y[^^,a^] represents the remaining terms in the Taylor expansion. By 
definition, these are at least cubic in ^ ^. Hence 



(63) 



exp(-r) = j P^-^exp {-S[a^] - L^[a^] • - ^L'^'K] • -V + T^'-^'^ 

We introduce a scaling parameter for the action, /i, (which we will set equal to one) via 

S ^ 

h 

r ^ -^r (64) 

We will show that an expansion in powers of h is consistent with the truncation used in 
deriving the moment hierarchy in section [3l The reason is that it organizes the terms in 
the expansion so that the true small parameters in the system, namely inverse distance from 
criticality and inverse number of inputs, are manifested. We thus consider an asymptotic 
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expansion T = Tq + hFi + h'^T2 + • • •. If we set Fq = S we obtain 

exp (^-ir) =exp (-ls[a^]^ J V^^exp (^-l-Li^-la^] . q, ^q>, - ^V[^ ^,a^] ■ + Oih) 

(65) 

Computing the corrections involves taking expectation values of the operator e~^l^ as well 
as other operators with respect to the Gaussian functional with covariance {L~^)^'^, which 
can be expanded as an infinite series of Gaussian moments. Fortunately, we can describe the 
terms in this series graphically using Feynman diagrams. A result of this analysis is that we 
can arrange the corrections to the effective action according to the number of loops in the 
Feynman diagrams, the order in h being given by the number of loops. |Zinn- Justin, 2002] . 
S[a^] provides the no loop, or "tree" level computation. The lowest order correction that these 
terms can produce is 0(1), which would be an 0(h) correction to the effective action (i.e. 1 
loop correction). This is because the corrections will be given by moments of operators which 
go as 1/h under a Gaussian functional distribution whose covariance goes as h. The terms Fi 
and higher produce still higher order corrections. We discuss in the appendix the connection 
between the h expansion, which is an artificial parameter, and the effective small parameters 
in the system, (i.e. the inverse of the distance to criticality and the inverse of the extent of 
connectivity within the network, as addressed in Sec 13. ip . 

To lowest order we obtain F[a^] = 5'[a^] which implies from (j42p that the equations of 
motion to lowest order are given by 

= {dt + a) a{x,t) - f {w * [d{x,t)a{x,t) + a{x,t)]) 

oa[x, tj 

- J (fx" d{x", t)f^^^ {w * [a{x", t)a{x", t) + a(x", t)] ) w{x" - x)a{x, t) = 

5S[a^] 



5a{x, t) 



{-dt + a) d{x,t) 
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/ 



d'^x" d{x", {w * [d{x", t)a{x", t) + a{x", t)] ) w{x" - x) [d{x, t) + 1] = 



from which we obtain a = (because there is no initial condition "source" term for d) 
and the mean field Wilson-Cowan equation ([1]). We can go to higher order by perform- 
ing a loop expansion on the path integral in ()65p and this next correction was computed in 
Buice and Cowan, 2007 . 

The importance of this approach is that as we consider successive orders in the loop ex- 
pansion, the effective action closes the system automatically. If we could calculate r[a^] for 
our model of interest, then we would have the exact equation of motion for the true mean of 
the theory. In essence, we are trading a closure problem for an approximation problem. The 
advantage gained is that we do not have to worry about the contributions of higher moments 
explicitly and we can consider explicitly the criteria allowing us to truncate the expansion. If 
there is an explicit loop expansion parameter, this truncation is straightforward. If not, as in 
our case, we must explicitly assess the regimes in which any truncation is valid. Even in the 
case where a truncation fails, the loop expansion can provide guidance in terms of identifying 
classes of diagrams (i.e. terms in the expansion) that are relevant in appropriate limits, which 
could be summed. 

We can now generalize this procedure for equations of motion for an arbitrary number 
of cumulants by considering a generating functional for an arbitrary number of "composite 
operators" [Cornwall et al., 1974| . In the case of the first and second cumulants, a{x,t) and 
C{x,y,t), we define the composite cumulant generating functional 




(66) 
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We now perform a double Legendre transform to obtain the effective composite action 



+\ [a^(x, t)a^{x', t') + hC^^{x, t; x' , t')] ■ t; x', t') 



r[a^,C^u] = -W[J^',K^''] + J^'ix,t)-a^ix,t) (67) 
1 
2 

with conditions 



6J^'{x,t) 



a^(x,t) (68) 



6W[J>^, _ 1 ^^^^^^ ^ ^^^^^^^ ^. ^gg^ 



6Kt'''{x,t;x',t') 2 



and 



'^^"'(x'^y^ = J'^(x, t) + ^a,{x', t') ■ [i^^-(x, t; x', t') + K'^-'ix', t'; x, t)] (70) 
dr[a^,C^.] _ (71) 



SC^i^{x,t', x' ,t') 2 

The equations of motion are obtained by setting = and K'^'"'^ = in (I70p and (j71|) . We 
calculate the equations of motion to lowest order for this system in the appendix. The results 
are 



(9j + a) a{x, t) — f {w -k a{x, t)) 
d'^x'd'^x" f^^\x,t)w{x-x")w{x-x')Cu{x',t;x",t) = (72) 



and for Cf^u, we get: 

{dt + a) Ci^-i{x,t;xo,to) — J dfix' f'^^^ {w -k a{x,t))w{x — x')Ci-i{x\t;xo 
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= 5{x-XQ)5{t-to) (73) 

{-dt + a) C_i,i(3;,t;xo,to) - j cfix' f^^^ {w-ka{x',t)) w{x' - x)C_i^i(x', t; xq, to) 

= S{x-xo)6{t-to) (74) 

{dt + a) Cii{x,t;xo,to) — J d'^x' f^^^ {w -k a{x,t)) w{x — x')Cii{x' ,t;xo,to) 

- J d'^x' f^^^wira{x,t))w{x - x')a{x',t') + f^^'^{w*a{x',t))w{x'-x)a{x,t) C_i,i(x', t; xq, to) 

= (75) 

together with the conditions 

Cii{x,t; x' ,t') = Cii{x' ,t';x,t) (76) 

C_i,i(x, t;x',t') = C7i,_i(x', t';x,t) (77) 

The 2-point correlation designated as Cfj_u{x, t; x', t') generahzes the cumulant C(x, y, t) in (fTSj) 
to include both the unequal time 2-point correlation (Cn) and the linear response (Ci^_i). 
The equation for C(x,y,t) is obtained by taking the equation for limt_,tQ Cii(x, t; xq, to) + 
Cii(xo, to; X, t) which results in (I30p . Note that we have also produced an equation for the 
Green's function of the theory Ci^_i as well as its time reversed counterpart. Time reversal 
involves swapping the fields (p and ^p. In the time reversed theory, (f plays the role of activity. 
Time reversal does not give an equivalent theory since our Markov process is not time reversal 
invariant in general. In field theory language, this system of equations is known as the 2PI 
equations and we adopt this moniker for conveniencqj- We can then analogously derive nPI 
equations for any number of normal ordered cumulants. 



i"2PI" stands for 2 Particle Irreducible. The effective action T[a^] is the generating functional of IPI graphs, 
which means that the graphs which determine T[a^] cannot be disconnected by cutting any single line of the graph. 
Similarly, P[a^,C^i/] is 2PI is the sense that graphs contributing to it cannot be disconnected by cutting two lines. 
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With the moment hierarchy approach, in order to produce better approximations we are 
required to truncate further in the hierarchy. This can quickly produce unwieldy equations. 
The loop expansion provides an alternative in that corrections to the generalized equations 
can be produced with a diagrammatic expansion, namely the one which calculates T2[a^, C^i,], 
from which the corrections to the equations can be calculated. 



5 All-to-All Networks, Finite Size Effects, and Sim- 
ulations 



We consider the example of an all-to-all system, wherein each neuron connects to the entire 
network. Mean field theory should work well in this case because the number of post-synaptic 
neurons reduces the coupling of the fluctuations. In this case, the fluctuations reasonably 
reduce to corrections due to the finite size of the network, as we would expect. We take 
the weight function to be a constant, normalized by the number of neurons in the system 
Wij = wq/N for some wq. The generalized Wilson-Cowan equations become 




(78) 




k 



k 




(79) 
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We can simplify this further by taking the initial conditions 0^(0) = and Cij{0) = and 
assuming homogeneous external input /. This corresponds to initial conditions in the network 
determined by a Poisson distribution with mean qq. Cjj(O) = — aj(0)5jj would indicate starting 
with precisely Oi(0) active neurons at i at time t = 0. Then symmetry reduces the equations 
to 

dta{t) = -aa{t) + f {woa{t) + I{t)) 

+^f"{woa{t)+Iit))wlC{t) (80) 



j^C{t) = -2aCit) + 2f'{woa{t) + I{t))woC{t) 

+^2f'{woa{t)+I{t))woait) (81) 



where we have defined 



i 



Note that as ^ cxd the source term for C{t) vanishes, which implies that C{t) decouples 
from the equation for o(t), which then reduces to the standard Wilson-Cowan equation. The 



matrix Tij in this case is the function 



r[a'']=a-f'{woa^ + I{t))wo (84) 
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The steady state value of C(t) is given by 

(B5) 

The relative size of the fluctuations in steady state is determined by the product iVr[ao]. Large 
networks or networks distant from a bifurcation experience reduced correlations. 

We now examine the phase plane of this simplified system. For concreteness, consider the 
firing rate function /(s) to be 

f{s) = tanh(s)e(s) (86) 

where G is the Heaviside step function. At mean field level (i.e. consider C{t) to be zero) 
equilibria are determined by solutions of the equation 



aao = /(ao) (87) 

Figure [1] graphically displays the solutions of this equation. From the figure, we see that the 
equation exhibits a bifurcation as the value of a (the slope of the straight line in the figure) 
is decreased. The critical value for this bifurcation is a = 1.0. We will refer to the non zero 
stable equilibrium as the "activated" state. The generalized activity equations ()80p and (j8ip 
also exhibit a bifurcation. The phase plane is shown in Figure [2] for a = 0.5,0.9 and 1.0 with 
N = 100 and in Figure E] for = 10 with the same values of a. As expected, the steady state 
value of C is larger for A^ = 10. Note that the nullclines for C{t) display divergences associated 
with r approaching zero. In addition note that the location of the bifurcation is different. For 
A^ = 100 the bifurcation happens near a = 0.9 and with A^ = 10 the bifurcation happens for 
0.9 < a < 1.0. Because the generalized equations are a coupled system, it is possible that more 
interesting bifurcation structure may be manifested. In addition to the fixed points which exist 
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Figure 1: Graphical depiction of solutions to equation (1571) for various values of a. 

in mean field theory there is a new fixed point. Whereas the bifurcation at mean field level is 
a pitchfork bifurcation, that of the generalized equations is a saddle node bifurcation, with an 
unstable fixed point emerging as the a and C nullclines cross. There is always a fixed point at 
a = and C = because it is an absorbing state. 

Importantly, we see that we can alter the bifurcation structure by adding a forcing or source 
term to the correlation function C{t) equation, linearly shifting the C nullcline. This removes 
the stable fixed point for high a (the one associated with the activated state in mean field 
theory). Because of this, we see that we can disrupt the activated state by stimulating the 
system with an input such that the correlation is sourced more strongly than the mean. We 
can use this to "turn off' the activated state by synchronizing the network. These correlations 
drive the system to the absorbing state of the full model. To reverse this deactivation, we 
simply drive the system with Poisson noise (i.e. force the equation for a{t) but not C{t)). 
Compare this to the effect demonstrated in [Laing and Chow, 2001| in which synchronized 
activity associated with fast synapses led to the destabilization of activity which the Wilson 
Cowan equation predicts to be stable. For a saturating firing rate function (more generally 
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Figure 2: Phase planes for the all-to-all generalized equations with a) a = 0.5, h) a = 0.9, and c) 
a = 1.0. = 100. Solid (black) lines are a nullclines; dotted (blue) lines are C nullclines. 
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Figure 3: Phase planes for the all-to-all generalized equations with a = 0.5 on the left, a = 0.9 in 
the center, and a = 1.0 on the right. = 10. Solid (black) lines are a nuUclines; dotted (blue) 
lines are C nullclines. 
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a function such that f"{s) < in the appropriate region) increased correlations inhibit the 
mean activity aj(t). 

We now demonstrate the utihty of the generahzed activity equations (j27p and (|28p for 
describing the full Markov system ^ away from a bifurcation point. In order to simulate 
the Markov model we use a form of the Gillespie algorithm and take expectation values over 
many time evolutions of the system. We use the firing rate function (j86p . It is important 
to remember when comparing results with simulations that we use the F{n) whose normal 
ordered form is f{s). In this all-to-all case, we need only consider the correction arising from 
the quadratic portion of the firing rate function since the corrections will go as powers of the 
weight function, which in this all-to-all case means they carry factors of In particular, 

to lowest order we have 



We plot ai{t) and Cii{t) versus t for various values of a and N in Figures [3] through [71 (Note 
that we are numerically evaluating the generalized equations, not the simplified equations in 
(|80p and (|8ip . and comparing them with Monte Carlo simulations of the Markov system; the 
plots overlay data for each of the N neurons.) We initialize the network with Poisson initial 
conditions: ai(0) = 2 and Cjj(O) = 0. The simulations of the Markov process are averaged 
over 10^ instances. One can see that the equations match the simulations quite well away 
from the critical point. As one approaches the bifurcation, however, the simulations begin 
to deviate. At a = 0.5, mean field and the generalized equations each match the simulated 
Markov process. As one approaches the mean field bifurcation point of a = 0.9, the mean 
field equations no longer match well, but the generalized equations account for the deviation. 
From Figure [5l one can see that as we approach a = 1.0, the estimate of the correlation from 
the generalized equations becomes poorer. 
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The plots for = 10 in Figures [6] and [7] demonstrate the breakdown of the generahzed 
equations. There is aheady a significant deviation of both mean field and the generalized 
equations at a = 0.5. Naturally, the discrepancy is accounted for by the poor estimation of 
the correlation at this level. As we near a = 1.0, the estimate of the correlations begins to 
grow, whereas the simulated correlation is dropping to zero. 

One can see that, even though the theory begins to deviate from the simulations near 
criticality, we still capture the loss of stability of the active state, even for N = 10. This is 
due to the growth of correlations which negatively feedback on the mean due to the negative 
second derivative of /(s). We can use this feature to observe the effect of correlated input 
directly by adding a term to the Markov process which provides multiplicative Gaussian noise. 
In particular, we add a transition at rate: 



where is a Gaussian noise source. The purpose of the step function is to prevent an individual 
neuron from getting a kick which will drive the activity negative. Note that we are not using 
an "input" term as we have defined it. Because the firing rate function is strictly positive, we 
cannot use a stimulus such that the mean transition rate is strictly zero. However, we can use 
a stimulus such that the source to the correlation function is much stronger than the mean. 
For our example, to maximize the effect, we do not use an input to the firing rate function but 
instead add another transition to the Markov process separate from the firing rate function and 
the decay transition. This allows us to source only the correlation function. Although this is 
artificial, this transition adds terms to the coupled equations which would approximate those 
from some input with zero mean and large correlations. In the following simulations we used 
a = 0.5, a = 100 and N = 100. The correlated input was given between t = 20 and t = 30. 




(89) 
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Figure 4: a{t) vs. t for a) a = 0.5, h) a = 0.9, and c) a = 1.0. = 100. Dotted (green) lines are 
solutions of mean field theory. Dashed (red) lines are solutions of the generalized equations. Solid 
(black) lines are expectations values of data from simulations of the Markov process. 
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Figure 5: C{t) vs. t for a) a = 0.5, h) a = 0.9, and c) a = 1.0. = 100. Dashed (red) lines 
are solutions of the generalized equations. Solid (black) lines are expectations values of data from 
simulations of the Markov process. 



44 



0.5 
0, 



0.5- 



'0 



10 



20 
t 



30 



40 



(a) 




(b) 



2 
1.5 




'0 10 20 30 40 



Figure 6: a{t) vs. t for a) a = 0.5, h) a = 0.9, and c) a = 1.0. = 10. Dotted (green) lines are 
solutions of mean field theory. Dashed (red) lines are solutions of the generalized equations. Solid 
(black) lines are expectations values of data from simulations of the Markov process. 
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Figure 7: C{t) vs. t for a) a = 0.5, b) a = 0.9, and c) a = 1.0. = 10. Dashed (red) lines 
are solutions of the generalized equations. Solid (black) lines are expectations values of data from 
simulations of the Markov process. 
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(a) ' (b) 

Figure 8: Response of all-to-all network to correlated input, a = 0.5, = 100. a) the response to 
correlated input with a = 100. b) the response to a Poisson process with rate A = 10. Note the 
change in scale between the two plots. 

In the absence of external correlated input, these parameters result in the active state being 
stable, as shown in Figure [H As can be seen in Figure [8l the use of the global noise source 
results in the "switch off' behavior predicted by the generalized equations. If we instead use 
a Poisson process to provide this external stimulation, one can also see in Figure H] that the 
network responds and then reverts back to the active state. The reason for the explicit shut 
off is that the system has an absorbing state to which it is driven. The more important point 
is that the correlated input is acting as a source of inhibition whereas the Poisson input serves 
as an excitatory input. A linear system or one in which the firing rate function /(s) is such 
that f"{s) > will not exhibit this behavior. Wc chose this particular example for sourcing 
Cij{t) in order to separate the effects of sourcing ai{t) as well. Given a more complicated noise 
source, one would need to examine the phase plane or solve the equations, after determining 
the effects of the noise source on the normal ordered cumulants. 
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6 Discussion 



We have demonstrated a formalism for constructing a minimal extension of a rate equation 
to include correlated activity. We have explicitly computed the lowest order results of this 
extension which provide coupled equations for the mean activity, two-point correlations, and 
linear responses. These results indicate that correlations can have an important impact on the 
dynamics of a rate equation, affecting both stability and the structure of bifurcations. Our 
argument relied upon inferring a "minimal" Markov process. Our use of the Doi-Peliti path 
integral formalism guides our assertion that our inferred Markov process is the simplest one 
compatible with both the rate equations and their interpretation as measuring some stochastic 
counting process. Thus, a general extension for any type of rate equation should share the 
same basic structure that we have described here. We performed this construction on a Markov 
process consistent with the Wilson-Cowan equation but our prescription would work equally 
well with any Markov process. 

In keeping with this idea, our results have something in common with other approaches to 
understanding correlations in neural networks. El Boustani and Destexhe |E1 Boustani and Destexhe, 2009] 
attempt to derive a Markov model for the asynchronous irregular states of an underlying neural 
system and explore the moment hierarchy of that Markov model. We take the opposite ap- 
proach, beginning with a presumed set of rate equations and asking what possible restrictions 
can be placed upon the correlation functions knowing only the dynamics of the rate model. 
Their hierarchy is truncated via scaling and finite size, whereas our hierarchy's truncation 
(and the truncation of the loop expansion) arises through the distance to a bifurcation in the 
rate equations. Ginzburg and Sompolinsky [Ginzburg and Sompolinsky, 1994] propose a sim- 
ple Markov model and study its moment hierarchy. For the correlations, they achieve results 
similar to the tree level of our loop expansion. An important point of departure is that we 
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consider the recurrent effects the correlations have upon the mean field, which we demonstrate 
can be sufficiently significant to alter the structure of the bifurcation. 

As we predict, our theory breaks down sufficiently close to a bifurcation. Examining the 
dynamics near the critical point requires a different form of analysis such as a renormalization 
argument. An example was presented in [Buice and Cowan, 2007| where it was argued that 
a large class of neural models would exhibit scaling laws near a bifurcation like those of 
the Directed Percolation model [Janssen and Tauber, 2005] . The predictions of this scaling 
coincide with measurements made in cortical slices of "avalanches" [Beggs and Plenz, 2003| . 



If criticality is important for neural function Beggs, 2008 , then these sort of approaches will 
be more important for future work and our rate model extension will be less applicable. 

In contrast, supporting the potential usefulness of our rate model extension is the fact that 
large neural connectivity suppresses correlations and aids the truncation of the hierarchy, an 
analogous result to the Ginsburg criterion in equilibrium statistical mechanics. In addition, we 
demonstrated that Poisson-like input in general pushes the system into configurations in which 
the correlations are suppressed relative to the mean. All of this suggests that our extension 
will be at least as applicable as the rate models themselves. 

Regarding that applicability, both the Markov process and the rate equations assume a 
large degree of underlying asynchrony in the network. The expansion we describe should 
be appropriate for networks in which a relatively small amount of synchrony at the level of 
individual neurons is developing. The coupled correlation function captures this effect. If the 
population is being dominated by neuron level synchrony, then the Markov process should no 
longer hold as a description of the system. Population level synchrony as captured by the 
original rate model, however, should have no effect on our analysis. In other words, there will 
be correlation effects acting on oscillating states, for example, such as presumably correlation 
induced modulation of the frequency of the oscillation. We will demonstrate this in future 
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work. 

An important outstanding point is that we have posited this Markov process based on 
the original interpretation of the Wilson-Cowan equations as dynamical equations for the 
probabilistic activity of neurons. Although our Markov process is the most "natural" given 
the transitions in the Wilson-Cowan equations, there is no a priori reason to suppose that 
this Markov process reflects the probabilistic dynamics of a physiologically based neural model 
or of real neurons precisely because there is nothing which mandates this interpretation. Per 
the renormalization analysis of [Buice and Cowan, 2007| , measuring scaling laws in cortex will 
provide a means of identifying classes of models (by identifying the relevant universality class) 
but this will in no way distinguish between models within the same class. Distinguishing 
models within the same class will require the measurement of non-universal quantities. This 
would likely mean some relatively precision measurements of response functions in cortical 
activity. 

Nonetheless, we feel our approach is a useful starting point for understanding effects be- 
yond mean field. We have demonstrated a correlation induced loss of stability in an all-to-all 
network. This effect should carry over to non-homogeneous solutions such as bump solutions 
or traveling waves. Likewise, correlations will modify important aspects of mean field solutions 
such as dispersion relations. Our approach enables this dispersion relation to be calculated. In 
addition to stability, our equations are a useful starting point towards understanding the wan- 
dering of bump solutions. They also provide a natural means of studying beyond mean field 
effects of correlation based learning. A model of visual hallucinations in cortex based on the 
Turing mechanism has explained many hallucinatory effects (such as the various Kluver form 
constants). Since the Turing mechanism is based upon bifurcations, it is an interesting question 
to what extent the coupling with correlations effects the results of the hallucination analysis. 
Our approach may provide this model with a means of explaining further hallucinations not 
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covered by the model in Bressloff et al., 2002 



It remains an important question how to connect our Markov and generahzed rate model 
approach with models of deterministic neurons. While the formalism admits almost any gain 
function, there remains the question of connecting this gain function to, for example, the 
transfer function for some neural model of which the Markov process is some approximation. 
This is of course not a question of the analysis of Markov models but of the applicability of 
rate models as high level descriptions of more detailed neural models. Answering this question 
will likely involve a kinetic theory formulation of the neural systems, such as the one pursued 
in [Hildebrand et al., 2006 Buice and Chow, 2007 . Having derived the generalized equations, 
it is also now important to explore their further consequences for phenomena such as pattern 
formation. There are also many avenues to extend this model and this approach. The Markov 
process can be enlarged to account for synaptic adaptation by adding a synaptic time variable 
to the neural configuration. Likewise noise in the transitions themselves, whether spatial or 
temporal, is easily incorporated into the action. A reduction of the resulting theory would no 
longer satisfy the Markov property, although there may be certain natural assumption (such 
as slow dynamics for the auxiliary field) that could allow one to regain Markovicity with an 
approximate model. This would allow us to construct extended Wilson-Cowan equations which 
incorporate these and other aspects of neural dynamics. These questions will be explored in 
future work. 
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A Composite Operator Effective Action and the 2PI 
equations 

Here we derive the 2PI equations. We begin with the generating functional: 



P^^exp ( -- 



^ j d'^x d'^x' dtdt' <^^{x,t)Kf"'{x,t;x',t')<^^{x',t') 



(90) 



where we have introduced a parameter h = 1 for bookkeeping purposes. The generaHzed 
effective action is given by 



-VF[J^,ii:n + j d'^xdtJ''{x,t)a^{x,t) (91) 
+ d'^xd'^x' dtdt' [af,{x,t)a^{x',t') + /iC^^(x, t; x', t')] K'"'{x,t;x' ,t') 



with and K'^'^ given by 



6W[J^',K''^] 
~~6J^^x^Tj~ 

6Ki"'{x,t;x',t') 



a^{x,t) 



1 



- [af,{x, t)au{x', t') + hCf,u{x, t; x' , t')] 



(92) 
(93) 



The path integral representation is thus 



exp(--r[a^,C^i,]) 



exp ( -- 



S[^f,]- J d'^xdt Jf'{x,t){<^f,{x,t)-af,{x,t)) 

^ j d'^xd'^x' dtdt' {<^>f,{x,t)^uix',t') - a^{x,t)a^{x' ,t') 
hC^,{x, t; x', i'))i^^'(x, t; x', t')] ) (94) 
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which we transform to a new variable = ^^{x,t) — a^{x,t), set 



Sr[a^^C^ + i / d^'x' dt' a,{x\ t') \K^'''{x, t; x' , t') + t' ■ x, t)\ (95) 

oa^ 2 J 



and expand S'f^^ + a^] = S[a^] + J d'^xdt{L^'[af,]^^ + {l/2)^^L^"'[af,]^^) + V['^^„a^] to obtain 

exp(-ir[a^,C7^,]) = exp{-^S[a^]) j P^^exp (^-i j d'^xdt (^L^^J^-^ + ^L^'^fa^]^^*, 
+ V[^^,a^]- J d''xdtTf''^[a^,C^,]^^{x,t) 

- \j d'^xd'^x' dtdt' {'iff,{x,t)^^{x',t') - hCf,^{x,t;x',t'))K'"'{x,t;x',t') 



where we use 



We consider the expansion F = Fq + hVi + /i^F2, where F2 contains all terms of order h? and 
higher. Setting Tq = S gives 



exp(-iFK,C^,]) = exp(--^5[o^]) j P^^exp (^-i j d'^x dt]^L>"'[a^,]^ , 
- j d''xdt{hT^^%^,C^,] + h^Tf[a^,C^,])^^{x,t) 

-\j d^x dt dt' {^^{x, t') - hC^,,{x, t; x',t'))Kf"'{x, t; x' , t') + V[^^, a^] 



We now fix K^'^ according to 



(5F[a^,C^i,j 



5Cfiv{x, t] x' , t') 



r°'^''[a^,C^,.] = -hK>^''{x,t;x',t') 



(96) 
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which gives 



exp(-irK,C^,]) = exp [-ls[af,] - ^Tr rO-^'^ia^, 



/ 



-i J d'^x d^x' dt dt' ^,{x, m,{x', tO^ ^^^^g ^^^^ + V[^,, a,] 



(97) 



where 



TiT^^^"'[a^„C^y\C^y = j d'^xd'^x'dtdf 



^Cui/{x, t; x' , t') 



-a 



(98) 



We need to extract the order h contributions from the functional integral. We expect the effect 
of r'^'^'' is to replace L^*^ with the full inverse two point function (C"^)^*^. This in turn will 
affect the normalization of the integral. Because of this we expect the order h contribution to 
the effective action to be: 



TiK, C^^] = ^Tr In {C~y + ^TrL^'''[a^]C^, + constant 



(99) 



Substituting this into expression ([U7|) gives us 



exp(-ir[a^, C^,]) = exp (-^5^] - ^TTL^'-[a^]C,,, + ^TrV " hTrT^^'^'ia^, C^,]C^, 



V^^ exp 



+y[*^,a^] - j d''xdt{hTf[a^,C^,] + h^Tf[a^,C^,])^^{x,t) 



-^j d'^xd'^x dtdt' ^^{x,t)^^{x',t')h 



(5C, 



(100) 



54 



We can extract the normalization of the functional integral using ilndetC^^ = -ilndet(C~i) 
and the identity In det ^ = Tr In A to obtain 



exp(-ir[a^,C^,]) = exp - - ^Trln {C 



detC^,. / P^-^exp 



+V[^^,a^] - J d^xdt{hrf[a^,C^,] + h^rf[a^,C^,])^^{x,t) 



(fx d'^x dt dt' 'i^f,{x,t)'^^{x',t')h- 



6C, 



The factor of the determinant serves as a normalization for the functional integral, which is 
now in a form that will only contribute to the effective action at order /i^. Thus we have 



rK, C^u] = S[a^] + ^/iTr In {0-^ + ^/iTrL^'^ia^JC^, 
+r2[a^, Cf„y] - ^/iTr Ifj,^ 



(102) 



Now we can calculate the equations of motion to a given loop order from equations ([70 
and (|7ip . Using equation (170p . the equations of motion for the mean field are 



5S[a 



(5r2[Q^,c^^j 



6a^ 2 6a^ 



6a,. 







(103) 



The equations for C^'^ are 



(104) 
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which we can invert to get 

L'^'ia.n. + f = V (105) 

In particular, if we ignore loop corrections (i.e. only consider first order in h, recalhng that r2 
is 0{h?')), we get 

C^, = L-l[a,,] (106) 

where L~l[a^] is the inverse of L^'^[a^]. In the absence of interactions, L~l is the two-point 
function, as expected. 

We can use the loop expansion to draw some conclusions about the applicability of per- 
turbation theory. Since T2[a^,C^u\ is second order and is the sum of vacuum two particle 
irreducible graphs, every graph contributing to it must be at least of two loop order. Every 
internal line represents a factor of C^j^u and so each graph contributing to r2 [a^ , Cfj_u\ niust 
have at least two factors of C^y, each of which will either be equal to 0, or be attenuated (in 
steady state) by the same exponents which attenuate the magnitude of C{x,y,t) away from 
a bifurcation, according to equations (jl20m22]) . Thus the argument that C^u is small away 
from the critical point extends to every term in the expansion for the generalized equations. 

The caveat here is that there is a class of diagrams which couple the lowest order expression 
for a given moment to the mean field. Although these graphs are suppressed by the distance 
to criticality, each of these is of the same order. We are assisted by two facts. The first 
is that the source terms for each of these moments at lowest order will be proportional to 
derivatives of the firing rate function. If f{s) is sufficiently smooth, this will suppress higher 
order contributions. In addition, each coupling will go as an additional factor of where 
Nm was defined in section [3TT] as the smallest number of inputs to any given neuron. Thus the 
connectivity in cortex will serve to "average out" sources to the mean from higher moments. 
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This will be the case as long as we can bound the total input to any given neuron. 



B Tree level equations of motion 

In order to calculate the expansion for the equations of motion, we need to compute the value 
of both L^'^ and We compute the lowest order correction here. 

First we find the intermediate results (which give us the classical equations of motion for 
a and a): 



5S[afj] 



{dt + a) a{x, t) — f {w -k [d{x, t)a{x, t) + a{x, t)]) 



^ 5d{x,t) 

d'^x" d{x", t)f^^^ {w ★ [d{x", t)a{x", t) + a{x", t)] ) w{x" - x)a{x, t) 



{-dt + a) d{x,t) 



^ 6a{x,t) 

d'^x" d{x", {w * [d{x", t)a{x", t) + a{x", t)] ) w{x" - x) [d{x, t) + 1] 



(107) 
(108) 



from which follows 



L-^'-^[af,]ix,t;x',t') 



L-^'^[a^]{x,t;x',t') 



- /(^) (x, t)w{x - x')a{x', t') + /(^) {x', t')w{x' - x)a{x, t)j 6{t - t') 

- j d{x", t)/(2) {x", t)w{x" - x)a{x, t)w{x" - x')a{x' , t)5{t - t') (109) 
{dt + a) 5{x - x')5{t - t') - f^^\x, t)w{ X — x') [a(x', t') + l] 8{t — t') 

d'^x"d{x",t)f^^\x",t")w{x" - x)6{x - x')5{t - t') 

d'^x"d{x\t)f^'^\x\t")w{x" - x)a{x,t)w{x" - x') [dix',t') + l] 6{t - t') 

(110) 

i^'"^ K] {x, t- x', t') = {-dt + a) 5{x - x')5{t - t') - /(^) {x', t)w{x' - x) [d{x, t) + 1] 5{t - t') 
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d'^x"a{x",t)f^^\x",t")wix" - x)6ix - x')5{t - t') 
d'^x"a{x",t)f^^\x",t")w{x" - x')a{x',t)w{x" - x) [d{x,t) + 1] 6{t - t') 

(111) 

d'^x"a{x",t)f'^^\x",t)w{x" - x)w{x" - x') [d{x,t) + 1] [d{x',t) + l] 6{t - t') 

(112) 



The terms f^'^\x,t) indicate the nth derivative of /. Note that we have suppressed the 
argument, so that /*^")(x,t) = Z*^"-* {w -k [d{x,t)a{x,t) + a{x,t)]) 

We can now write down the equations of motion from (jlOSp . minus the loop corrections. 
The first "diagonal" equation (for (—1,-1)) is: 



(at + a)Ci,_i(x,t;xo,to) - / d'^x' f^^\x,t)w{x - x') [d{x',t') + l] Ci,_i(x', t; xq, to) 

d'^x" d{x", t)/(^) {x", t")w{x" - x)Ci^-i{x, t- xo, to) 
d'^x" d'^x' d{x'\t)f^'^\x" ,t")w{x" - x)a{x,t)w{x" - x') [d{x\t') + l] Ci,_i(x', t; xq, to) 



d'^x' 



f^^\x,t)w{x - x')a{x',t') + f^^\x',t')w{x' - x)a(x,t)j C_i _i(x', t; xq, to) 



d'^x" d'^x' a(x", t)/(2) {x", t)w{x" - x)a(x, t)w{x" - x')a{x\ t)C„i,_i (x', t; xq, to) 

= 5{x - xo)5{t - to) 



(113) 



The second "diagonal" equation (for 11): 



{-dt + a) C_i,i(x, t; XQ, to) - / d'^x' f^^\x', t)w{x' - x) [a(x, t) + 1] C_i,i(x', t; xq, to) 

d'^x" d{x", t)/(^) (x", t")w{x" - x)C7_i,i(x, t; xo, to) 
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■ I d'^x" <fx' o(x", (x", t")w{x" - x')a{x\ t)w{x" - x) [a{x, t) + 1] C_i,i(x', t; xq, to) 
d-^x' a(2;", t)/^^) {x" , t)w{x" - x)w{x" - x') [d{x, t) + 1] [d{x' , t) + l] Cn {x' , t; xq, to) 



5(x - xo)6{t - to) 



(114) 



The "off-diagonal" equations are (starting with —1,1): 



{dt + a)Cu{x,t;xo,to)- / d'^x' f^^\x,t)w{x - x') [d{x' ,t') + l] Cu{x' ,t;xo,to) (115) 

d'^x" d{x", t)/(^) (x", t")''^(^" - x)Cii{x, t; xo, to) 
dVdVo(x",t)/(2)(x",t")w;(x" -x)a(x,t)w;(x" - x') [o(x',t') + l] Cii(x', t; xq, to) 



d'^x' 



f^^\x,t)w{x - x')a{x',t') + /(^)(x',t')w(x' - x)a(x,t) C_i,i(x', t; xq, to) 
dVd'^x'a(x",t)/(2)(x",t)u;(x" -x)a(x,t)u;(x" -x')a(x',t)C7_i,i(x',t;xo,to) = 

and the other (1, —1): 

(-9t + a)C7_i,_i(x,t;xo,to) - / (iV /^(x', t)u;(x' - x) [a(x, t) + 1] C_i,_i(x', t; xo, to) (116) 

d'^x" a(x", t)/(^) (x", t")w{x" - x)C_i,_i (x, t; xo, to) 
dVdVa(x",t)/(2)(x",t")w^(x" -x')a(x',t)w;(x" -x) [a(x,t) + 1] C„i,_i(x', t; xo, to) 



The "mean field" portion of the equations of motion (jl03p are obtained from equations (jl07p 
and (jl09p (by setting the LHS to zero) . The remainder of the equations of motion are "classi- 
cal" terms dependent on the correlation functions, and loop corrections. The latter are given 
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by the term ^h-g^TrL^" 



[o-iilC^u in equation (jl03|) . The term in the trace is, of com'se, the 



sum of the LHS of equations 11131 and 11141 

We can simphfy the equations for the mean field by reahzing that any term involving C_i^i 
or Ci-i can be ignored because they will only appear in the form C_i^i(x', t; x, t), i.e. at equal 
initial and final times. These will be zero. This can be seen as either the "initial condition" 
for the linear response terms or as a manifestation of the "e(0)" problem in quantum field 
theory. [Zinn-Justin, 2002] 

Furthermore, we can use some results from the full theory. In particular, we have 



because causality enforces that (f operators can't contract with anything "in the past". 
The equation for a{x, t) is then: 



C_i _i(x,t;x',t') = 



(117) 



a{x,t) = 



(118) 



{dt + a) a{x, t) — f {w -k a{x, t)) 



J 



cfix' d'^x" f^'^\x,t)w{. 



x-x")w{x-x')Cii{x',t;x",t) = 



(119) 



Applying these same simplifications to the equations for C^y, we get: 





x - x')Ci _i(x',t;xo,to) 



5{x - xo)6{t - to) 



(120) 





6{x - xo)S{t - to) 



(121) 





x - x')Cii{x',t;xo,to) 
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f^^^ {w -k a{x,t)) w{x — x')a{x' ,t') + f^^^[w-ka{x',t))w{x' — x)a{x,t) C_i^i(x', t; j;o> *o) 

= (122) 



together with the conditions 



Cii{x,t;x',t') = Cii{x',t']x,t) 
C-i^i{x,t;x',t') = Ci-i{x',t';x,t) 



(123) 
(124) 



References 

[Abbott and van Vreeswijk, 1993] Abbott, L. and van Vreeswijk, C. (1993). Asynchronous 
states in networks of pulse-coupled oscillators. Phys. Rev. E, 48(2): 1483-1490. 

[Amari, 1975] Amari, S. (1975). Homogeneous nets of neuron-like elements. Biological Cyber- 
netics, 17:211-220. 

[Amari, 1977] Amari, S. (1977). Dynamics of pattern formation in lateral-inhibition type 
neural fields. Biological Cybernetics, 27:77-87. 

[Amit and Brunei, 1997a] Amit, D. and Brunei, N. (1997a). Model of global spontaneous 
activity and local structured activity during delay periods in the cerebral cortex. Cereb. 
Cortex, 7(3):237-252. 

[Amit and Brunei, 1997b] Amit, D. J. and Brunei, N. (1997b). Dynamics of a recurrent net- 
work of spiking neurons before and following learning. Network: Computation in Neural 
Systems, 8(4):373-404. 

[Amit et al., 1985] Amit, D. J., Gutfreund, H., and Sompolinsky, H. (1985). Spin-glass models 
of neural networks. Physical Review A, 32(2): 1007-1018. 

61 



[Beggs, 2008] Beggs, J. M. (2008). The criticality hypothesis: how local cortical networks 
might optimize information processing. Philosophical Transactions of The Royal Society A, 
366(1864) :329-343. 

[Beggs and Plenz, 2003] Beggs, J. M. and Plenz, D. (2003). Neuronal avalanches in neocortical 
circuits. The Journal of Neuroscience, 23(35):11167-11177. 

[Beshel et al., 2007] Beshel, J., Kopell, N., and Kay, L. (2007). Olfactory bulb gamma oscil- 
lations are enhanced with task demands. Journal of Neuroscience, 27(31):8358-65. 

[Bressloff et al., 2002] Bressloff, P., Cowan, J. D., Golubitsky, M., Thomas, P. J., and Wiener, 
M. C. (2002). Geometric visual hallucinations, euclidean symmetry and the functional 
architecture of striate cortex. Philosophical Transactions of the Royal Society of London B, 
356:299-330. 

[Brunei and Hakim, 1999] Brunei, N. and Hakim, V. (1999). Fast global oscillations in net- 
works of integrate-and-fire neurons with low firing rates. Neural Computation, 11(7):1621- 
1671. 

[Brunei and Hansel, 2006] Brunei, N. and Hansel, D. (2006). How noise affects the synchro- 
nization properties of recurrent networks of inhibitory neurons. Neural Comput., 18(5):1066- 
1110. 

[Buice and Chow, 2007] Buice, M. A. and Chow, C. C. (2007). Correlations, fluctuations, 
and stability of a finite-size network of coupled oscillators. Physical Review E (Statistical, 
Nonlinear, and Soft Matter Physics), 76(3):031118. 

[Buice and Cowan, 2007] Buice, M. A. and Cowan, J. D. (2007). Field theoretic approach to 
fluctuation effects for neural networks. Physical Review E, 75:051919. 

[Cai et al., 2004] Cai, D., Tao, L., Shelley, M., and McLaughlin, D. W. (2004). An effective 
kinetic representation of fluctuation-driven neuronal networks with application to simple 

62 



and complex cells in visual cortex. Proceedings of the National Academy of Sciences of the 
United States of America, 101(20):7757-7762. 

[Cohen and Grossberg, 1983] Cohen, M. and Grossberg, S. (1983). Absolute stability and 
global pattern formation and parallel storage by competitive neural networks. IEEE Trans- 
actions on Systems, Man, and Cybernectics, 13:815-826. 

[Coombes, 2005] Coombes, S. (2005). Waves, bumps, and patterns in neural field theories. 
Biological Cybernetics, 93:91-108. 

[Cornwall et al., 1974] Cornwall, J. M., Jackiw, R., and Tomboulis, E. (1974). Effective action 
for composite operators. Physical Review D, 10(8):2428-2445. 

[Doi, 1976a] Doi, M. (1976a). Second quantization representation for classical many-particle 
system. Journal of Physics A: Mathematical and General, 9(9):1465-77. 

[Doi, 1976b] Doi, M. (1976b). Stochastic theory of diffusion controlled reaction. Journal of 
Physics A: Mathematical and General, 9(9): 1479-95. 

[El Boustani and Destexhe, 2009] El Boustani, S. and Destexhe, A. (2009). A master equa- 
tion formalism for macroscopic modeling of asynchronous irregular activity states. Neural 
Computation, 21:46-100. 

[Ermentrout, 1998] Ermentrout, B. (1998). Neural networks as spatio-temporal pattern- 
forming systems. Reports on Progress in Physics, 61(4):353-430. 

[Ermentrout and Cowan, 1979] Ermentrout, G. B. and Cowan, J. D. (1979). A mathematical 
theory of visual hallucination patterns. Biological Cybernetics, 34:137-150. 

[Fourcaud and Brunei, 2002] Fourcaud, N. and Brunei, N. (2002). Dynamics of the firing 
probability of noisy integrate-and-fire neurons. Neural Comput., 14(9):2057-2110. 



63 



[Fusi and Mattia, 1999] Fusi, S. and Mattia, M. (1999). Collective behavior of networks with 
linear (vlsi) integrate-and-fire neurons. Neural Computation, ll(3):633-652. 

[Gerstner, 1995] Gerstner, W. (1995). Time structure of the activity in neural network models. 
Phys. Rev. E, 51(l):738-758. 

[Gerstner, 2000] Gerstner, W. (2000). Population dynamics of spiking neurons: Fast tran- 
sients, asynchronous states, and locking. Neural Computation, 12:43-89. 

[Gerstner and Kistler, 2002] Gerstner, W. and Kistler, W. M. (2002). Spiking Neuron Models. 
Cambridge University Press. 

[Ginzburg and Sompolinsky, 1994] Ginzburg, I. and Sompolinsky, H. (1994). Theory of corre- 
lations in stochastic neural networks. Physical Review E, 50(4):3171-3191. 

[Golomb and Hansel, 2000] Golomb, D. and Hansel, D. (2000). The number of synaptic inputs 
and the synchrony of large, sparse neuronal networks. Neural Comput., 12(5):1095-1139. 

[Gray and Singer, 1989] Gray, C. M. and Singer, W. (1989). Stimulus-specific neuronal oscil- 
lations in orientation columns of cat visual cortex. Proceedings of the National Academy of 
Sciences of the United States of America, 86(5):1698-1702. 

[Hildebrand et al., 2006] Hildebrand, E. J., Buice, M. A., and Chow, C. C. (2006). A kinetic 
theory of coupled oscillators. arXiv, nlin.CD. 

[Hopfield, 1982] Hopfield, J. (1982). Neural networks and physical systems with emergent 
collective computational abilities. Proceedings of the National Academy of Sciences, 79. 

[Hopfield, 1984] Hopfield, J. (1984). Neurons with graded response have collective computa- 
tional properties like those of two state neurons. Proceedings of the National Academy of 
Sciences, 81. 



64 



[Ichimaru, 1973] Ichimaru, S. (1973). Basic principles of Plasma Physics: A Statistical Ap- 
proach. W.A. Benjamin Advanced Book Program. 

[Janssen and Tauber, 2005] Janssen, H.-K. and Tauber, U. C. (2005). The field theory ap- 
proach to percolation processes. Annals of Physics, 315(1):147 - 192. Special Issue. 

[Laing and Chow, 2001] Laing, C. and Chow, C. (2001). Stationary bumps in networks of 
spiking neurons. Neural Computation, 13(7): 1473-1494. 

[Ly and Tranchina, 2007] Ly, C. and Tranchina, D. (2007). Critical analysis of dimension 
reduction by a moment closure method in a population density approach to neural network 
modeling. Neural computation, 19(8):2032-92. 

[Martin et al., 1973] Martin, P. C, Siggia, E. D., and Rose, H. A. (1973). Statistical dynamics 
of classical systems. Physical Review A, 8(l):423-437. 

[Mattia and Del Giudice, 2002] Mattia, M. and Del Giudice, P. (2002). Population dynamics 
of interacting spiking neurons. Phys. Rev. E, 66(5):051917. 

[N., 2000] N., B. (6 May 2000). Dynamics of sparsely connected networks of excitatory and 
inhibitory spiking neurons. Journal of Computational Neuroscience, 8:183-208(26). 

[Nicholson, 1992] Nicholson, D. R. (1992). Introduction to Plasma Theory. Krieger Publishing 
Company. 

[Nykamp and Tranchina, 2000] Nykamp, D. Q. and Tranchina, D. (2000). A population den- 
sity approach that facilitates large-scale modeling of neural networks: Analysis and an 
application to orientation tuning. J. Comp. Neurosci, 8. 

[Peliti, 1985] Peliti, L. (1985). Path integral approach to birth-death processes on a lattice. 
Journal de Physique, 46:1469-1483. 



65 



[Peretto, 1984] Peretto, P. (1984). Collective properties of neural networks: a statistical 
physics approach. Biological Cybernetics, 50:51-62. 

[Plesser and Gerstner, 2000] Plesser, H. E. and Gerstner, W. E. (2000). Noise in integrate- 
and-fire neurons: From stochastic input to escape rates. Neural Comput., 12(2):367-384. 

[Salinas and Sejnowski, 2002] Salinas, E. and Sejnowski, T. J. (2002). Integrate-and-fire neu- 
rons driven by correlated stochastic input. Neural Comput., 14(9):2111-2155. 

[Softky and Koch, 1993] Softky, W. R. and Koch, C. (1993). The highly irregular firing of 
cortical cells is inconsistent with temporal integration of random epsps. The Journal of 
Neuroscience, 13(l):334-350. 

[Soula et al., 2006] Soula, H., Beslon, G., and Mazet, O. (2006). Spontaneous dynamics of 
asymmetric random recurrent spiking neural networks. Neural Comput., 18(l):60-79. 

[Soula and Chow, 2007] Soula, H. and Chow, C. C. (2007). Stochastic dynamics of a finite-size 
spiking neural network. Neural Comput., 19(12):3262-3292. 

[Tauber et al., 2005] Tauber, U. C, Howard, M., and Vollmayr-Lee, B. P. (2005). Applica- 
tions of field-theoretic renormalization group methods to reaction&:ndash;diffusion problems. 
Journal of Physics A: Mathematical and General, 38(17):R79-R131. 

[Treves, 1993] Treves, A. (August 1993). Mean-field analysis of neuronal spike dynamics. 
Network: Computation in Neural Systems, 4:259-284(26). 

[Van Vreeswijk and Sompolinsky, 1996] Van Vreeswijk, C. and Sompolinsky, H. (1996). Chaos 
in neuronal networks with balanced excitatory and inhibitory activity. Science, 274:1724- 
1726. 

[Wilson and Cowan, 1972] Wilson, H. R. and Cowan, J. D. (1972). Excitatory and inhibitory 
interactions in localized populations of model neurons. Biophysical Journal, 12. 

66 



[Wilson and Cowan, 1973] Wilson, H. R. and Cowan, J. D. (1973). A mathematical theory of 
the functional dynamics of cortical and thalamic nervous tissue. Kybernetik, 13:55-80. 

[Zinn- Justin, 2002] Zinn- Justin, J. (2002). Quantum Field Theory and Critical Phenomena. 
Oxford Science Publications, 4 edition. 



67 



